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Abstract 



We present the program SecDec 2.0 which contains various new features: First, 
i-S^ . it aUows the numerical evaluation of multi-loop integrals with no restriction on 

the kinematics. Dimensionally regulated ultraviolet and infrared singularities are 
isolated via sector decomposition, while threshold singularities are handled by a 
(— I , deformation of the integration contour in the complex plane. As an application we 

present numerical results for various massive two-loop four-point diagrams. SecDec 
2.0 also contains new useful features for the calculation of more general parameter 

■ integrals, related e.g. to phase space integrals. 
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Nature of problem: 

Extraction of ultraviolet and infrared singularities from parametric integrals ap- 
pearing in higher order perturbative calculations in gauge theories. Numerical inte- 
gration in the presence of integrable singularities (e.g. kinematic thresholds). 

Solution method: 

Algebraic extraction of singularities in dimensional regularisation using iterated 
sector decomposition. This leads to a Laurent series in the dimensional regularisa- 
tion parameter e, where the coefficients are finite integrals over the unit-hypercube. 
Those integrals are evaluated numerically by Monte Carlo integration. The inte- 
grable singularities are handled by choosing a suitable integration contour in the 
complex plane, in an automated way. 

Restrictions: Depending on the complexity of the problem, limited by memory and 
CPU time. The restriction that multi-scale integrals could only be evaluated at Eu- 
clidean points is superseded in version 2.0. 

Running time: 

Between a few minutes and several days, depending on the complexity of the prob- 
lem. 
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1 Introduction 



Currently we are in tlie fortunate situation of being confronted with a wealth of 
high energy collider physics data, enabling us to test our present understanding 
of fundamental interactions and to explore physics at the TeV scale. However, 
the accuracy which has been or will be reached by the experiments has to be 
matched by comparable precision in the theory predictions, and in most cases 
this means that calculations beyond the leading order in perturbation theory 
are necessary. 

It is well known that in the calculation of higher order corrections, various 
types of singularities can arise at intermediate stages of the calculation. For 
example, loop integrals can contain ultraviolet (UV) as well as infrared (IR) 
singularites, phase space integrals over unresolved massless particles lead to in- 
frared singularities, and there can be integrable singularities due to kinematic 
thresholds. The UV and IR singularities can be regularised by dimensional 
regularisation, such that they appear as poles in 1/e, which cancel when the 
different parts of the calculation are combined to a physical observable. How- 
ever, before such cancellations are possible, the 1/e poles have to be extracted. 
In the calculation of multi-loop integrals or real radiation at higher orders, 
this usually leads to the task of factorising the poles from complicated multi- 
parameter integrals. Sector decomposition [T][2|3] is a method to achieve such 
a factorisation. The program SecDec 1.0, presented in performs this task 
in an automated way. Other public implementations of sector decomposition 
can be found in [SfGfTfH] . see also [9]. The method already has been applied 
in various calculations, listing all of them is beyond the scope of this paper, 
for a review see [10]. Here we just mention that there are also fruitful combi- 
nations of sector decomposition with other techniques, e.g. non-linear trans- 
formations [TT], Mellin-Barnes and differential equation techniques [7]|12|13] . 
high-energy expansions P^fT^] , or in the context of subtraction for unresolved 
double real radiation at NNLO [TM71[TM9]|20]|2T]|22]|23]I^ A method devel- 
oped over many years |25f26f27f28f29f30] to calculate one- and two-loop inte- 
grals numerically in the physical region also partly uses sector decomposition, 
in combination with a careful analysis of the singularity structure of certain 
classes of integrals and the use of functional relations between loop integrands. 

A limitation of the program SecDec 1.0 was the fact that the numerical 
integration of multi-scale integrals was only possible for Euclidean points, or, 
more precisely, values of the Mandelstam invariants and masses for which 
the denominator of the integrand is guaranteed to be of definite sign. For 
physical applications which go beyond one-scale problems, it is however crucial 
to be able to deal with integrable singularities, usually related to kinematic 
thresholds, in addition to the singularities in e. The program SecDec 2.0 
is able to achieve this task, by an automated deformation of the integration 
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contour into the complex plane. This allows the numerical calculation of multi- 
scale integrals in the physical region in an automated way. Non-planarity of 
the considered integral does not add any extra complications. Adding more 
mass scales also does not necessarily increase the complexity of the calculation 
with this method, as additional masses usually lead to a simpler IR singularity 
structure and therefore lead to less functions in the iterated decomposition. 
Therefore one can for instance obtain numerical results for two-loop integrals 
involving several mass scales where analytic methods reach their limit. 

The method of contour deformation in a multi-dimensional parameter space 
in the context of perturbative calculations has been pioneered in [31] and 
later has been refined in various ways to be applied to calculations at one 
loop |32f33f34f35f36f37f38f39] and at two loops [iO|41f42|H3] . 
Another purely numerical method uses an extrapolation from large to small 
values of the (analytically infinitesimal) parameter moving the integration con- 
tour away from poles on the real axis P^|l5] . Numerical methods using dis- 
persion relations, differential equations and/or numerical integration of Mellin- 
Barnes representations also have been worked out, see e.g. |46|47|48|49|50|51|52j 
However, most numerical methods to calculate multi-scale integrals beyond 
one loop so far are either limited to specific types of integrals, or the parame- 
ters for the numerical integration have been carefully adapted to the individual 
integrals by the authors. 

The aim of the work presented here is to provide a public program where the 
user can calculate multi-scale integrals without worrying too much about the 
details of the integrand. The singularity structure does not have to be known 
beforehand (but certainly the user has to make sure that, after the extraction 
of the poles regulated by dimensional regularisation, only integrable singular- 
ities remain). The program contains a sophisticated procedure to check and 
adjust the contour deformation parameters to optimize the convergence. For 
complicated integrals, the convergence can nonetheless depend critically on 
the settings for the numerical integration; therefore we also offer the possibil- 
ity for the user to choose various parameters at the input level to tune the 
deformation. 

We should note that, even though the functions which are produced after the 
factorisation of the singularites in e are available in algebraic form, they are 
usually too complicated to be integrated analytically. Therefore the final in- 
tegration is done by the Monte Carlo methods, meaning that the precision 
which can be achieved is limited, but in favour of a gain in general applica- 
bility. We also should remark that the method is applicable to any number of 
loops in principle, however memory problems in the algebraic part where the 
functions are generated, or bad numerical convergence can be expected if the 
complexity of the integral is very high. 

The structure of this paper is as follows. In Section [21 we briefly describe 
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the general framework. Section [3] gives an overview of the structure of the 
program and the new features introduced in SecDec version 2.0. Section 
m contains installation and usage instructions, while examples and results 
are presented in Section El A brief user manual is given in the Appendix. 
Detailed documentation is also coming with the code which is available at 
http : / / secdec . hepf orge . org. 



2 General framework 



The procedure of factorising endpoint singularities from parameter integrals by 
iterated sector decomposition is described in Here our main concern is the 
numerical integration for physical kinematics after the endpoint singularities 
have been extracted. Our method to do so is based on contour deformation, 
described in detail in Section [2^ The following section serves to introduce 
some basic concepts. 



2.1 Feynman integrals 



We choose a scalar integral for ease of notation. Tensor integrals only lead 
to an additional function of the Feynman parameters and invariants in the 
numerator. For more details we refer to [IfTU] . 

A scalar Feynman integral G m. D dimensions at L loops with propagators, 
where the propagators can have arbitrary, not necessarily integer powers z/j, 
has the following representation in momentum space: 



G=/nd%^, — ^ 

d% = ^d%, P,{{k},{p},m]) = q]-m] + i5, (1) 
in 2 

where the qj are linear combinations of external momenta Pi and loop momenta 
ki. Introducing Feynman parameters leads to 
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J^{x) = det (M) Qj^ji^Qi- J -'^^ 



(2) 
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W(x)=det(M), = Y^J ■ 



The functions U and J-" also can be constructed from the topology of the 
corresponding Feynman graph [53|54fl0] , and the implementation of this con- 
struction in SecDec2.0 is one of the new features of the program. 

For a diagram with massless propagators, none of the Feynman parameters 
occurs quadratically in the function J-" = J-'q . If massive internal lines are 

present, J-" gets an additional term J^{x) = J^o{x) + U{x) Xjiri]. 

W is a positive semi-definite function. A vanishing U function is related to 
the UV subdivergences of the graph. Overall UV divergences, if present, will 
always be contained in the prefactor T{Njj—m — LD/2). In the region where all 
invariants formed from external momenta are negative, which we will call the 
Euclidean region in the following, J-" is also a positive semi-definite function 
of the Feynman parameters Xj. Its vanishing does not necessarily lead to an 
IR singularity. Only if some of the invariants are zero, for example if some 
of the external momenta are light-like, the vanishing of J-" may induce an IR 
divergence. Thus it depends on the kinematics and not only on the topology 
(like in the UV case) whether a zero of J-" leads to a divergence or not. The 
necessary (but not sufficient) conditions for a divergence are given by the 
Landau equations [55|56ll57j : 




If all kinematic invariants formed by external momenta are negative, the nec- 
essary condition J-" = for an IR divergence can only be fulfilled if some of 





(3) 



(4) 
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the parameters Xi go to zero. These endpoint singularities can be regulated 
by dimensional regularisation and factored out of the function J-" using sector 
decomposition. The same holds for dimensionally regulated UV singularities 
contained in U. However, after the UV and IR singularities have been ex- 
tracted as poles in 1/e, for non- Euclidean kinematics we are still faced with 
integrable singularities related to kinematic thresholds. How we deal with these 
singularities will be described in the following section. 



2.2 Deformation of the integration contour 



Im(z) 




Re(z) 



Fig. 1. Schematic picture of the closed contour avoiding poles on the real axis. 

Unless the function J-" in eq. ([2]) is of definite sign for all possible values of 
invariants and Feynman parameters, the denominator of a multi-loop integral 
will vanish within the integration region on a hypersurface given by the so- 
lutions of the Landau equations. If eq. ([3]) has a solution where Xj > OVj, 
all particles in the loop go simultaneously on-shell. This corresponds to a 
leading Landau singularity, which is not integrable (for real values of masses 
and momenta). However, the integrand can also diverge for certain values of 
kinematical invariants and Feynman parameters which represent a subleading 
Landau singularity, corresponding to an integrable singularity of logarithmic 
or squareroot type, related to normal thresholds. In these cases, we can make 
use of Cauchy's theorem to avoid the non-physical poles on the real axis by a 
deformation of the integration contour into the complex plane. As long as the 
deformation is in accordance with the causal i6 prescription of the Feynman 
propagators, and no poles are crossed while changing the integration path, the 
integration contour can be changed such that the convergence of the numerical 
integration is assured. Using the fact that the integral over the closed contour 
in Fig. [T]is zero, we have 
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X{z{x)) , 



(5) 



where the Xj are real, while Zi are complex, describing a path parametrized by 
the variables Xi. The i5 prescription for the Feynman propagators tells us that 
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the contour deformation into the complex plane should be such that the imag- 
inary part of T should always be negative. For real masses and Mandelstam 
invariants Sjj, the following Ansatz pT|33f34] is therefore convenient: 

z{x) = X — i r (x) 

Tk = Xxk{l - Xk) ^ . (6) 

dxk 

The derivative of J-" in eq. ([6]) is smallest in the extrema and largest where the 
slope is maximal. Hence, unless we are faced with a leading Landau singularity 
where both J-" and its derivatives with respect to Xi vanish, the deformation 
leads to a well behaved integral at the points where the function J-" vanishes. 
A closed integration contour is guaranteed by the factors Xk and (1 — Xk), 
keeping the endpoints fixed. In terms of the new variables, we thus obtain 

T{z{x)) = J-(f) - z A E x,{l - X,) + 0{\') , (7) 

such that J-" acquires a negative imaginary part of order A. Hence, the size of 
A determines the scale of the deformation. More technical details about the 
deformation are given in Section 13.31 



2.3 Parameter integrals 



The program SecDec can also factorise singularities from parameter integrals 
which are more general than the ones related to multi-loop integrals. The only 
restrictions are that the integration domain should be the unit hypercube, and 
the singularities should be only endpoint singularities, i.e. should be located at 
zero or one. Contour deformation is not available in the subdirectory general, 
because the sign of the imaginary part telling us how to deform the contour is 
not fixed a priori for general functions, in contrast to loop integrals. However, 
we plan to implement the use of Cauchy's theorem where applicable in a 
future version. Currently we assume that the singularities are regulated by 
non-integer powers of the integration parameters, where the non-integer part 
is the e of dimensional regularisation or some other regulator. The general 
form of the integrals is 

1 ™ 
dXi... dXN\{P^{x,{a]T , (8) 

where Pi{x,{a}) are polynomial functions of the parameters Xj, which can 
also contain some symbolic constants {a}. The user can leave the parameters 
{a} symbolic during the decomposition, specifying numerical values only for 
the numerical integration step. This way the decomposition and subtraction 




8 



steps do not have to be redone if the values for the constants are changed. 
The Vi are powers of the form z/j = + h^e (with such that the integral is 
convergent). Note that half integer powers are also possible. 



3 The SecDec program 



3. 1 Structure 



The program consists of two parts, an algebraic part and a numerical part. 
The algebraic part uses code written in Mathematica [SS| and does the de- 
composition into sectors, the subtraction of the singularities, the expansion 
in e and the generation of the files necessary for the numerical integration. 
In the numerical part, Fortran or C++ functions forming the coefficient of 
each term in the Laurent series in e are integrated using the Monte Carlo in- 
tegration programs contained in the CuBA library [59p0] . or BASES [M]. The 
different subtasks are handled by perl scripts. The flowchart of the program 
is shown in Fig. [2] for the basic building blocks to calculate multi-loop inte- 
grals. To calculate parameter integrals which are not necessarily related to 
loop integrals, the structure is the same except that contour deformation is 
not available. For more details about the features in the part general we refer 
tog]. 
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Fig. 2. Flowchart showing the main steps the program performs to produce the 
result Laurent series in e. 

The directories loop and general have the same global structure, only some 
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of the individual files are specific to loop integrals or to more general para- 
metric functions. The directories contain a number of perl scripts steering the 
decomposition and the numerical integration. The scripts use perl modules 
contained in the subdirectory perlsrc. 

The Mathematica source files are located in the subdirectories src/deco (files 
used for the decomposition), src/subexp (files used for the pole subtraction 
and expansion in e) and src/util (miscellaneous useful functions). The doc- 
umentation, created by robodoc [02] is contained in the subdirectory doc. It 
contains an index to look up documentation of the source code in html format 
by loading masterindex.html into a browser. 

In order to use the program, the user only has to edit the following two files: 

• par am. input: (text file) 

specification of paths, type of integrand, order in e, output format, param- 
eters for numerical integration, further options 

• Template .m: (Mathematica syntax) 

■ for loop integrals: specification of loop momenta and propagators, resp. 
of the topology; optionally numerator, non-standard propagator powers, 
space-time dimensions 

■ for general functions: specification of integration variables, integrand, vari- 
ables to be split 

The program comes with example input and template files in the subdirectories 
loop/demos respectively general/demos, described in detail in 

3.2 New features of the program 

Version 2.0 of SecDec contains the following new features, which will be 
described in detail in this section, while examples will be given in Section [51 

• loop part: 

■ Multi-scale loop integrals can be evaluated without restricting the kine- 
matics to the Euclidean region. This has been achieved by performing a 
(numerical) contour integration in the complex plane. The program auto- 
matically tries to find an optimal deformation of the integration path. 

■ For scalar multi-loop integrals, the integrand can be constructed from the 
topological cuts of the diagram. The user only has to provide the vertices 
and the propagator masses, but does not have to provide the momentum 
flow. 

■ The files for the numerical integration of multi-scale loop diagrams with 
contour deformation are written in C-|--|- rather than Fortran. For inte- 
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grations in Euclidean space, both the Fortran and the C++-versions are 
supported. The choice between Fortran and C++ can be made by the 
user in the p ar am . input file by choosing either language=Cpp (default) 
or language=fortran . 

■ A parallelisation of the algebraic part for Mathematica versions 7 and 
higher is possible if several cores are available. 

■ The most recent version of the Cuba library, CuBA-3.0 (beta), is added 
to the program and used by default. The older version CuBA-2.1 is still 
supported. 

■ The rescaling of the kinematic invariants is now possible by choosing 
rescale=l (default is 0). 

• general part: 

■ The user can define additional (finite) functions at a symbolic level and 
specify them only later after the integrand has been transformed into a 
set of finite parameter integrals for each order in e. 

• both parts: 

■ The possibility to loop over ranges of parameter values is automated. 
In the following we describe the new features in more detail. 



3.3 Multi-scale loop integrals: Implementation of the contour deformation 

As explained in Section 12.21 singularities on the real axis can be avoided by 
a deformation of the integration contour into the complex plane. The overall 
size of the deformation is controlled by the parameter A defined in eq. (jH]). 

The convergence of the numerical integration can be improved significantly 
by choosing an "optimal" value for A. Values of A which are too small lead to 
contours which are too close to the poles on the real axis and therefore lead to 
bad convergence. Too large values of A can modify the real part of the function 
to an unacceptable extent and could even change the sign of the imaginary 
part if the terms of order A^ get larger than the terms linear in A. This would 
lead to a wrong result. Therefore we implemented a four-step procedure to 
optimize the value of A, consisting of 

• ratio check: To make sure that the terms of order A^ in eq. ([7]) do not spoil 
the sign of the imaginary part, we evaluate the ratio of the terms linear and 
cubic in A for a quasi-randomly chosen set of sample points to determine 
the maximal allowed A = Xmax- 

• modulus check: The imaginary part is vital at the points where the real part 
of J-" is vanishing. In these regions, the deformation should be large enough 
to avoid large numerical fluctuations due to a highly peaked integrand. 
Therefore we check the modulus of each subsector function J^i at a number 
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of sample points, and pick the fraction of the value of Xmax which maximises 
the minimum of the modulus of J-i, i.e. the value of lambda which keeps J-i 
furthest from zero. 

• individual X{i,j) adjustments: If the values of |^ are very different in mag- 
nitude, it can be convenient to have an individual parameter X{i,j) for each 
subsector function J-'j and each Feynman parameter Xj. 

• sign check: After the above adjustments to A have been made, the sign of 
lm(J-') is again checked for a number of sample points. If the sign is ever 
positive, this value of A is disallowed. 

The contour deformation can be switched on or off by choosing 
contourdef= True/False in the input file par am. input. Obviously, the calcula- 
tion takes longer if contour deformation is done, so if the integrand is known 
to be positive definite, contour deformation should be switched off. We also 
should emphasize that for integrands with a complicated threshold structure, 
the success of the numerical integration can critically depend on the param- 
eters which tune the deformation, and on the settings for the Monte Carlo 
integration. In order to allow the user to tune the deformation, the following 
parameters can be adjusted by the user in the input file: 

lambda: the program takes the A value given in param . input as a starting 
point. If, after the program has performed the checks listed above, this A 
is found to be unsuitable or suboptimal, the value of A will be changed 
automatically by the program. The default is lamhda=1.0. 

largedefs: If the integrand is expected to have (integrable) endpoint singu- 
larities at Xj = or 1, the deformation should be large in order to move 
the contour away from the problematic region. If largedefs=l (default is 0), 
the program tries to enlarge the deformation at the endpoints. Further, the 
program performs a remapping to separate endpoint problems for xj — )■ 1 
from the ones for xj — )■ if this flag is set to one. This increases the number 
of functions, but it can be very useful to improve convergence. 

smalldefs: If the integrand is expected to be oscillatory and hence sensitive to 
small changes in the deformation parameter A, choosing the flag smalldefs=l 
(default is 0) will minimize the argument of each subsector function J-^ by 
varying A(i, j). 

3.4 Topology-based construction of the integrand 

As already mentioned in section HI the functions U and J-' can be constructed 
from the topology of the corresponding Feynman graph [53|54fl0] . without 
the need to assign the momenta for each propagator explicitly. The user only 
has to label the external momenta and the vertices. If an external momentum 
Pi is part of a vertex, this vertex needs to carry the label i. The labelling of 
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vertices containing only internal lines is arbitrary. In Template. m, the user 
has to specify proplist as a list of entries of the form {mslk], {i, j}}, where 
ms[fc] is the mass squared of the propagator connecting vertex i and vertex 
j. The mass label k must correspond the the kth entry of the list of masses 
given in param. input. While k needs to be the number labelling the masses, 
ms[/c] (with k being an integer) can be left symbolic during the decomposition. 
However, if the mass is zero, one has to put {0, because this changes 

the singularity structure at decomposition level. 

An example is given below, more examples can be found in the mathemat- 
ica template files template? 126 .m, templateBnp6* .m, template JapNP .m, 
templateggtt* .m in the subdirectory loop/demos. This feature of construct- 
ing the graph topologically is only implemented for scalar integrals so far. 
The original form of specifying the propagators by their momenta, as done in 
SecDecI.O, is still operational. The topology based construction is selected 
by defining cutconstruct= 1 in the input file. 



3. 5 Looping over ranges of parameters 

As the algebraic part can deal with symbolic expressions for the kinematic 
invariants or other parameters contained in the integrand, the decomposition 
and subtraction parts only need to be done once for the calculation of many 
different numerical points. Therefore it is desirable to automate the calculation 
of many numerical points to minimize the effort for the user. This is done 
using the perl script multinumerics.pl. The user should create a text file 
multiparamf ile in myworkingdir, and specify a number of options: 

• paramfile=myparamfile: specify the name of the parameter file. 

• pointname=myprefix: points calculated will have the names myprefixl, mypre- 
fix2, . . . 

• lines: the number of points you wish to calculate - if omitted all points 
(listed in separate lines) will be calculated. 

• xplot the number of the column containing the values which should be used 
on the X-axis of the plot (default is 1). 

After these options, the numerical values of the parameters for each point to 
calculate should be specified. In the loop directory, the number of values given 
for Sij, p1 and needs to be specified by numsij=, numpi2= and numms2=. 
An example can be found in loop/demos/multiparam. input. The following 
example explains how the numerical values for each point are written down in 
the general directory. If you wished to calculate three numerical points for 
a function where the symbols a, h (defined as symbols in the parameter input 
file) should take on the values (a, 6) = (0.1, 0.1), (0.2, -0.4), (-0.3, 0.9) then 
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the inputs in multiparamf ile for this would be: 

0.1,0.1 

0.2,-0.4 

-0.3,0.9 

Furthermore, one may wish to calculate the integrand for values of parameters 
at incremental steps. This is allowed, and the syntax is as follows: Suppose you 
wish to calculate each combination of s = 0.1, 0.2, 0.3 and t = 0.1, 0.3, 0.5, 0.7. 
The input for this is 

1711^013 = 0.1,0.1 

maxvals=0.3,0. 7 
stepvals=0.1,0.2 

Non-equidistant step values are also possible. For instance, to calculate every 
combination of a = 0.1, 0.2, 0.4, h = 0.1, 0.3, 0.6 the syntax would be: 
valuesl=0. 1,0.2,0.4 
values2=0. 1,0.3,0.6 

Please note that valuesl must appear before values2 in multiparamf ile. 
Examples can be found in general/demos/mult iparam. input or 
loop/demos/mult iparam. input. In the loop directory, there is a perl script 
helpmulti.pl which can be used to generate the files mult iparam. input 
automatically to avoid typing large sets of numerical values. 
In order to execute the script multinumerics .pi, the Mathematica-generated 
functions must already be in place. The simplest way to do this is to run the 
launch script, with exeflag=l in your parameter file. Then issue the command 
'./multinumerics. pi [-d myworkingdir -p multiparamfile] '. In single-machine 
mode {dusterflag=0) all integrations will then be performed, and the results 
collated and output as files in the directory specified in myparamf ile. In batch 
mode you will need to run the script again, with the argument '1', to collect 
the results, i.e. './multinumerics. pi 1 [-d myworkingdir -p multiparamfile] '. 
The script generates a parameter file for each numerical point calculated. To 
remove these intermediate parameter files (your original myparamfile will not 
be removed), issue the command './multinumerics. pi 2 [-d myworkingdir -p 
multiparamfile] This should only be done after the results have been collated. 

3.6 Leaving functions implicit during the algebraic part 

This feature is available in the part general to evaluate general parametric 
functions, where it is possible to include a "dummy" function depending on 
(some of) the integration parameters, the actual form of the function being 
specified only later at the numerical integration stage. There are a number of 
reasons why one might want to leave functions implicit during the algebraic 
stage. For example, squared matrix elements typically contain large but finite 
functions of the phase space variables in the numerator, so the algebraic part 
of the calculation will be quicker and produce much smaller intermediate files 
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if these functions are left implicit. Also, one might like to use a number of 
measurement functions and be able to specify or change them without having 
to redo the decomposition. To use this option, the Mathematica template file 
can contain a function which is left undefined, but needs to be listed under 
the option dummys in the parameter input file. Note that one may use more 
than one implicit function at a time, and that these functions can have any 
number of arguments. If symbolic parameters are also used, these do not need 
to be arguments of the implicit function. 

Once the template and parameter files are set up, the functions need to be 
defined explicitly so that they can be used in the calculation. The simplest way 
to do this is to prepare a Mathematica syntax file for each implicit function 
specified, and place them in the outputdir specified in your parameter file. 
Suppose you have a function named duml, a function of two variables, defined 
as duml{xi,X2) = 1 + Xi+X2- Then you should create a file duml . m, and insert 
the lines: 

intvars = {zi, Z2}] 
duml = 1 + Zi + Z2', 

where Zi, Z2 can be replaced by any variable name you wish, as long as they are 
used consistently in duml . m. Notice that for every function specified in dummys 
in your parameter file, there must be a Mathematica file dummyname.m with 
the correct name and syntax in the results directory. Once these Mathematica 
files are in place, issue the command 

^ createdummyfortran.pl [-d myworkingdir -p myparamfilef from the general 
directory. This generates the fortran files for the functions you defined, which 
are found in the same subdirectory as the originals. 

Of course you might prefer to write these fortran files yourself instead of 
having them generated by the program. This is certainly possible, however 
we recommend that you use this perl script to generate functions with the 
necessary declarations and then edit these. 

An example of this can be found in general/demos, with the files 
paramdummy . input , templatedummy .m, and the directory /testdummy. 



4 Installation and usage 

4.1 Installation 

The program can be downloaded from 
http : / / secdec . hepf orge . org. 

Unpacking the tar archive via tar xzvf SecDec.tar.gz will create a directory 
called SecDec with the subdirectories as described above. Then change to the 
SecDec directory and run ./install. 
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Prerequisites are Mathematica, version 6 or above, perl (installed by default 
on most Unix/Linux systems), a Fortran compiler (e.g. gfortran, ifort) or a 
C++ compiler if the C++ option is used. 



4-2 Usage 

(1) Change to the subdirectory loop or general, depending on whether you 
would like to calculate a loop integral or a more general parameter inte- 
gral. 

(2) Copy the files param. input and template .m to create your own param- 
eter and template files myparamfile . input, mytemplatef ile .m. 

(3) Set the desired parameters in myparamfile . input and specify the inte- 
grand in mytemplatef ile .m. 

(4) Execute the command ./launch -p myparamfile. input -t mytemplatefile.m 
in the shell. 

If you omit the option -p myparamfile. input., the file param. input will be 
taken as default. Likewise, if you omit the option -t mytemplatefile.m., the 
file template. m will be taken as default. If your files myparamfile. input, 
mytemplatefile.m are in a different directory, say, myworkingdir, use the 
option -d myworkingdir, i.e. the full command then looks like ./launch 
-d myworkingdir -p myparamfile. input -t mytemplatefile.m, executed from 
the directory SecDec/loop or SecDec/general. 

(5) Collect the results. Depending on whether you have used a single machine 
or submitted the jobs to a cluster, the following actions will be performed: 

• If the calculations are done sequentially on a single machine, the results 
will be collected automatically (via results .pi called by launch). The 
output file will be displayed with your specified text editor. 

• If the jobs have been submitted to a cluster, when all jobs have fin- 
ished, execute the command ./results. pi [-d myworkingdir -p myparam- 
file]. This will create the files containing the final results in the graph 
subdirectory specified in the input file. 

(6) After the calculation and the collection of the results is completed, you 
can use the shell command ./launchclean [graph] to remove obsolete files. 



It should be mentioned that the code starts working first on the most compli- 
cated pole structure, which takes longest. This is because in case the jobs are 
sent to a cluster, it is advantageous to first submit the jobs which are expected 
to take longest. 
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5 Examples and Results 



5.1 Massive two-loop integrals 

5.1.1 A two-loop three-point function 

In this example, we will demonstrate three of the new features of the SecDec 
2.0 program: the construction of J-", W directly from the topology of the graph, 
the evaluation of the graph in the physical region, and how results for a whole 
set of different numerical values for the invariants can be produced and plotted 
in an automated way. We will use the two-loop diagram shown in Fig. [3] as 
an example. Numerical results for this diagram have been produced in [2Z1E3], 
and an analytical result can be found in [61], where the diagram is called Pi2%- 




Fig. 3. Two- loop vertex graph Pi2q, containing a massive triangle loop. Solid lines 
are massive, dashed lines are massless. The vertices are labeled to match the con- 
struction of the integrand from the topology as explained in the text. 

The template file templateP126 .m in the demos subdirectory contains the fol- 
lowing lines: 

proplist={{ms[l],{3,4}},{ms[l],{4,5}},{ms[l],{5,3}},{0,{l,2}},{0,{l,4}},{0,{2,5}}}; 
onshell={ssp[l]-)- 0,ssp[2]-)- 0,ssp[3]-> sp[l,2]}; 

where each entry in proplist corresponds to a propagator of the diagram; 
the first entry is the mass of the propagator, and the second entry contains 
the labels of the two vertices which the propagator connects. The labels for 
the vertices are as shown in Fig.|3l Note that if an external momentum p^ is 
flowing into the vertex, the vertex must also have the label k. For vertices 
containing only internal propagators the labeling is arbitrary. The on-shell 
conditions in the above example state that Pi= P2 = 0, p\ = Si2 = s. Results 
for the e° part of graph Pi2e are shown in Fig. HI 

To run this example, from the loop directory, issue the command . /launch -d 
demos -p paramP 1 26 . input -t templateP126 .m. The timings for the finite 
part and a relative accuracy of 1%, using CuBA-3.0 [60], are about 15.5 sees 
for a typical point far from the s = 4m^ threshold, and 20.5 sees for a point 
close to threshold (s/m^ = 3.9) on an Intel(R) Core(TM)2 Quad CPU Q9400 
at 2.66GHz. 
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Fig. 4. Comparison of analytic and numerical results for the diag ram P\2Q using 
= 1. 

Producing data files for sets of numerical values 



To loop over a set of numerical values for the invariants s and once the 
C++ files are created, issue the command 

perl multinumerics.pl -d demos -p mult iparamPl 26 . input. This will run 
the numerical integrations for the values of s and specified in the file 
demos/multiparamP126 . input. The files containing the results will be found 
in demos/21oop/P126, and the files p-2.gpdat, p-l.gpdat and pO.gpdat 
will contain the data files for each point, corresponding to the coefficients of 
e~^, and e° respectively. These files can be used to plot the results against 
the analytic results using gnuplot. This will produce the files P126R0.ps, 
P126I0.ps which will look like Fig. H 
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5.1.2 Non-planar massive two-loop four-point functions 
The graph B^^ 




Fig. 5. The non-planar 6-propagator graph B^^. 
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Fig. 6. Results for the finite part of the graph Bq^ (a) with p'f and p'2 off-shell, (b) 



with mi = m2 = = me = 0.25, ms = m^ = 0. The imaginary part of Bq ' is 
zero in the range shown above. 



Next, we consider the non-planar 6-propagator two- loop four-point diagram 
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shown in Fig.|5l For light-like legs and massless propagators, the analytic result 
has been calculated in [65], where the graph is called B^^ . Here we give results 
for this topology for the cases where 

(a) p\ and p\ are off-shell 

(b) nil = ^2 = ^5 = 7^ 0, m3 = = 0. 

It is interesting to note that B^^'°' with p\ and p| being off-shell contains 
poles starting from while for light-like legs the leading pole is only 
due to cancellations related to the high symmetry of the graph. For Bq^'^, 
i.e. the graph with mi = m2 = = ^ 0, the leading pole is 1/e. Results 
for the finite parts of Bq^'"" and B^^'^ are shown in Figs.E] and [71 As in 




Fig. 7. The non-planar 6-propagator graph B^^'^ as a function of Si2 in a region 
containing a threshold. 

[65], an overall prefactor of r(l + 2e)r(l - e)Vr(l - 3e)/(l + 4e) has been 
extracted. For Fig.|6]we have used the numerical values Si2 = 3 while scanning 
over S23- For all the values given, S13 is determined by the physical constraint 
S12 + Sis + -523 = Pi + pi- For Bq^'"" we have set pf = pi = 1, while for B^^'^ , 
mi = m2 = m^ = m^ = 0.25 has been used. Fig.[7] shows Bq^'^ as a function 
of S12 with S23 = —0.4 and mi = m2 = m^ = m^ = 0.25. The numerical 
accuracy is about one permil, therefore the error bars are barely seen in the 
figures. 



The graph J 

In this example we consider a 7-propagator non-planar two-loop box integral 
where all propagators are massive, using mi = m2 = m^ = mg = m, 777,3 = 
7774 = 7777 = M, Pi=P2=P3=p'i = The labelling is as shown in Fig. [HI 

Numerical results for this integral have been calculated in [15] using a method 
based on extrapolation in the i5 parameter. Our results for = 50, M = 
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Fig. 8. Labeling of the masses for the non-planar graph J^^ . 

90,523 = —10"^ are shown in Fig. [9] and agreement with ref. |15] has been 
verified. The timings for the longest subfunction (both real and imaginary 
part) with a relative accuracy of one permil vary between about 20 seconds 
for a point far from threshold and about 500 seconds close to threshhold. 




Fig. 9. Results for a non-planar 7-propagator graph J^^ with all propagators mas- 
sive using m = 50, M = 90 and t = —10^. 



The graph ggtti 

The graph shown in Fig.[TU]occurs in the calculation of the two-loop corrections 
to heavy quark production. Numerical results for the two-loop amplitude in 
the qq initiated channel have been calculated in [51]. Analytic results in the qq 
channel and for some colour structures in the gg channel have been calculated 
in [66|67|I68] . Numerical results for the amplitude in the gg channel in the 
approximation s ^ have been calculated in [69)170] . 

For the individual graph shown in Fig.[TOl numerical results at Euclidean 
points have been given in [1]. Here we give numerical results for the non- 
planar master integral ggtti in the non-Euclidean region. For the results shown 
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in Fig. [m we used S23 = —0.4, S13 = — 0.1,p| = pI = —0.17, = m2 = 0.17. 
The analytical result for this master integral is not known yet. 




Fig. 10. Non-planar graph ocurring in the calculation of gg — > tt at NNLO. Blue 
(solid) lines denote massive particles. 
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Fig. 11. Results for the non-planar graph ggtt\. 



5.2 Defining implicit functions 



This example demonstrates the ability to leave certain functions implicit until 
numerical integration. Suppose we want to integrate 

f{x) = (xi + X2)~'^~'^''x^^~^^dumi{xi, X2, X3, XiY~^''dum2{x2, X4)^~^''cut{x3) 

with 



dumi{xi, X2, X3, Xi) = 2 + xl + xl + xl + xl + 4x1X2X3X4, — xfx^x^x^ 
dum2{x2,X4) = x\ + x\ + 15"^ + 4x2X4 — ^Jx^x^^ + 3x1x1 , 

where /3 is a symbol defined in the parameter file, durrii and dum2 should 
always be functions which cannot increase the singular behaviour of the in- 
tegrand, and so quantitative knowledge of their exact form is not required to 
guide the decomposition. Thus they can be left implicit, and only introduced 
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at the numerical integration stage. The function cut{x) is aheady defined by 
default as cut{x) = G(x — cutval), where the value of cutval is given in the 
parameter file. The command ../launch -p paramdummy. input -t template- 
dummy. m from the folder general/demos runs this example. The fortran files 
containing the explicit form of the functions durrii, dum2, cut are found in 
demos/test dummy. Another simple example for a dummy function would be a 
jet function as used in the JADE algorithm in e~^e~ annihilation, as described 
e.g. in [7l]. Input files for such an example are params23s35 JADE, input, 
templates23s35JADE.m in the general/demos directory. 



6 Conclusions 

We have presented the program SecDec2.0, which can be used to factorise 
dimensionally regulated singularities and numerically calculate multi-loop in- 
tegrals in an automated way. As a new feature of the program, it now can deal 
with fully physical kinematics, i.e. is not restricted to the Euclidean region 
anymore. A new construction of the integrand, based entirely on topological 
rules, is also included. The new features are demonstrated by several examples, 
among them a massive two-loop four-point function which is not yet known 
analytically. In addition, the program can produce numerical results for more 
general parameter integrals, as they occur for example in phase space integrals 
for multi-particle production with several unresolved massless particles. The 
program also offers the possibility to include symbolic functions which can be 
used for instance to define measurement functions like jet algorithms in a fiex- 
ible way. The program setup is such that the evaluation of several functions in 
parallel can lead to a major speed-up. To calculate full two-loop amplitudes 
involving several mass scales, the timings still leave room for improvement, 
but considering the fact that the method is very suitable for intense paral- 
lelisation, we are convinced that the program will be a very useful tool for a 
multitude of applications to higher order corrections in quantum field theories. 
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A User manual 



We list here all possible input parameters for the parameter file * . input and 
the Mathematica input file * . m. These two files serve to define the integrand 
and the parameters for the numerical integration. We describe here the exam- 
ple of loop diagrams; the input files in the subdirectory general to compute 
more general parametric functions is very similar. 



A.l Program input parameters 

This input file should be called *. input. The following parameters can be 
specified 

subdir suhdir specifies the name of the subdirectory to which the graph 
should be written to. If not yet existent, it will be created. The specified 
suhdir contains the directory specified in outputdir. 

outputdir The name for the desired output directory can be given here. If 
outputdir is not specified, the default directory for the output will have the 
graph name (see below) appended to the directory suhdir, otherwise specify 
the full path for the Mathematica output files here. 

The output directory will contain all the files produced during the decom- 
position, subtraction, expansion and numerical integration, and the results. 
The output of the decomposition into sectors is found in the outputdir di- 
rectly. The functions from subtraction and expansion and the respective 
files for numerical integration are found in subdirectories. The latter are 
named with the pole structure and contain subdirectories named with the 
respective Laurent coefficient. 

graph The name of the diagram or parametric function to be computed is 
specified here. The graph name can contain underscores and numbers, but 
should not contain commas. 

propagators The number of propagators the diagram has is specified here 
(mandatory) . 

legs The number of external legs the diagram has is specified here (manda- 
tory). 

loops The number of loops the diagram has is specified here (mandatory). 

cutconstruct If the graph to be computed corresponds to a scalar integral, 
the integrand (F and U) can be constructed via topological cuts. In this 
case set cutconstruct=l, the default is =0. If cutconstruct is switched on, 
the input for the graph structure (*.m file) is just a list of labels connecting 
vertices, as explained in Section 13.41 and IA.2I 

epsord The order to which the Laurent series in e should be expanded, start- 
ing from g-»^'^2:po«e^ pg^j^ specified here. The default is epsord=0 where the 
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Laurent series is cut after finite part e°. If epsord is set to a negative value, 
only the pole coefficients up to this order will be computed, 
prefactorfiag Possible values for the prefactorflag are (default), 1 and 2. 

• 0: The default pref actor (— l)'^r[A^ — N loops * Dim/ 2] is factored out of 
the numerical result. 

• 1: The default prefactor (— 1)^ r[A^ — Nloops * Dim/ 2] is included in the 
numerical result. 

• 2: Give the desired prefactor to be factored out in prefactor. 
prefactor If option 2 has been chosen in the prefactorflag, write down the 

desired prefactor in Mathematica syntax here. In combination with options 
or 1 in the prefactorflag this entry will be ignored Use Nn, Nloops and 
Dim to denote the number of propagators, loops and dimension (4-2eps by 
default). 

IBPfiag Set IBPflag=0 if integration by parts should not be used, =1 if it 
should be used. IBPflag=2 is designed to use IBP when it is more efficient to 
do so, and not otherwise. Using the integrations by parts method takes more 
time in the subtraction and expansion step and generally results in more 
functions for numerical integration. However, it can be useful if (spurious) 
poles of type x~'^~'^^ are found in the decomposition, as it reduces the power 
of X in the denominator. 

compiler Set a Fortran compiler (tested with gfortran, ifort, g77) if lan- 
guage= Fortran. Left blank, the default is gfortran. 

exefiag The exeflag is set to decide at which stage the program terminates: 

• 0: The iterated sector decomposition is done and the scripts to do the 
subtraction, the expansion in epsilon, the creation of the Fortran/C++ 
files and to launch the numerical integration are created (scripts batch* 
in the subdirectory graph) but not run. This can be useful if a cluster is 
available to run each pole structure on a different node. 

• 1: In addition to the steps done in 0, the subtraction and epsilon expansion 
is performed and the resulting functions are written to Fortran/C++ files. 

• 2: In addition to the steps done in 1, all the files needed for the numerical 
integration are created. 

• 3: In addition to the steps done in 2, the compilation of the Fortran/C++ 
files is launched to make the executables. 

• 4: In addition to the steps done in 3, the executables are run, either by 
batch submission or locally. 

clusterfiag The clusterflag determines how jobs are submitted. Setting clus- 
terflag=0 (default) the jobs will run on a single machine, setting it =1 the 
jobs will run on a cluster (a batch system to submit jobs). 

batchsystem If a cluster is used {clusterflag =1), this fiag should be set to 
to use the setup for the PBS (Portable batch system). If the fiag is set to 1 a 
user-defined setup is activated. Currently this is the submission via condor, 
but the user can adapt this to his needs by editing perlsrc/makejob.pm. 

maxjobs When using a cluster, specify the maximum number of jobs allowed 
in the queue here. 
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maxcput Specify here the estimated maximal CPU time (in hours). This op- 
tion is used to send a job to a particular queue on a batch system, otherwise 
it is not important. 

pointname The name of the point to calculate is specified here. It should 
be either blank or a string and is useful to label the result files in case of 
different runs for different numerical values of the Mandelstam variables, 
masses etc. 

sij The values for Mandelstam invariants Sij = {pi + pj)"^ in numbers are 
specified here (mandatory). The Sij should be < in the Euclidean region. 

pi2 The off-shell legs pf, p\.,... are specified here (mandatory), pj should be 
< in the Euclidean region. 

ms2 Specify the masses of propagators m\, m\,... here using the notation 
ms[i] for m? (mandatory). The masses should not be complex numbers. 

integrator The program for numerical integration can be chosen here. BASES 
[integrator =0) can only be used in the Fortran version. Vegas {integra- 
tor=l), Suave [integrator = 2) , Divonne [integrator =3, default) and Cuhre 
[integrator =4) are part of the CuBA library and can be used in both the 
Fortran and the C-|--|- version. In practice, Divonne usually gives the fastest 
results when using the C+-f- version. In the following we therefore concen- 
trate on the adjustment of the parameters needed for numerical integration 
using Divonne. For more details about the Cuba parameters we refer to 

m- 

cubapath The path to the Cuba library can be specified here. The default 
directory is [your path to SecDec]/ Cuba- 3.0. CuBA-3.0 is the newest version 
of the Cuba library and uses parallel processing during the numerical eval- 
uation of the integral. The older version (Cuba-2.1) is still supported and 
can be used. 

mcixeval Separated by commas and starting with the lowest order coefficient 
in e, specify the maximal number of evaluations to be used by the numerical 
integrator for each order in e. If maxevalis not equal to mineval, the maximal 
number of evaluations does not have to be reached. 

mineval Separated by commas and starting with the lowest order coefficient 
in e, specify the number of evaluations which should at least be done before 
the numerical integrator returns a result. The default is 0. 

epsrel Separated by commas and starting with the lowest order coefficient in 
e, specify the desired relative accuracy for the numerical evaluation. 

epsabs Separated by commas and starting with the lowest order coefficient in 
e, specify the desired absolute accuracy for the numerical evaluation. This 
becomes useful in the cases where the integrated result is close to zero. 

cubafiags Set the cuba verbosity fiags. The default is 2 which means, the 
Cuba input parameters and other useful information, e.g. about numerical 
convergence, are echoed during numerical integration. 

keyl Separated by commas and starting with the lowest order coefficient in e, 
specify keyl which determines the sampling to be used for the partitioning 
phase in Divonne. With a positive keyl, a Korobov quasi-random sample of 
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keyl points is used. A keyl of about 1000 (default) usually is a good choice. 

key2 Separated by commas and starting with the lowest order coefficient 
in e, specify key2 which determines the sampling to be used for the final 
integration phase in Divonne. With a positive key2, a Korobov quasi-random 
sample is used. The default is key2=l which means, the number of points 
needed to reach the prescribed accuracy is estimated by Divonne. 

key3 Separated by commas and starting with the lowest order coefficient in 
e, specify the keyS to be used for the refinement phase in Divonne. Setting 
key3=l (default), each subregion is split once more. 

maxpass Separated by commas and starting with the lowest order coefficient 
in e, specify how good the convergence has to be during the partitioning 
phase until the program passes on to the main integration phase. A maxpass 
of 3 (default) is usually sufficient to get a quick and good result. 

border Separated by commas and starting with the lowest order coefficient 
in e, specify the border for the numerical integration. The points in the 
interval [0, border] and [1 — border, 1] are not included in the integration 
but are extrapolated from points further from the endpoints. This can be 
useful if the integrand is known to be peaked at endpoints of the integration 
variables. 

maxchisq Separated by commas and starting with the lowest order coeffi- 
cient in e, specify the maximally allowed the end of the numerical 
integration. 

mindeviation Separated by commas and starting with the lowest order co- 
efficient in e, specify the deviation two sample averages in one region can 
show without being treated any further. 

These parameters are advanced options 

primarysectors Specify a list of primary sectors to be treated here. If left 
blank, primarysectors defaults to all, i.e. 1 to the number of propagators, 
will be taken. This option is useful if a diagram has symmetries such that 
some primary sectors yield the same result. 

multiplicities Specify the multiplicities of the primary sectors listed above. 
List the multiplicities in same order as the corresponding sectors above. If 
left blank, default multiplicities (=one) are set automatically. 

infinitesectors A list of primary sectors to be redone differently because they 
lead to infinite recursion can be specified here, infinitesectors must be left 
empty for the default strategy to be applied. 

togetherfiag This flag deflnes whether to integrate subsets of functions for 
each pole order separately togetherflag=0{deia.u\t) or to sum all functions 
for a certain pole order and then integrate togetherflag=l. The latter will 
allow cancellations between different functions and thus give a more realistic 
error, but should not be used for complicated diagrams where the individual 
functions are large already. 

editor Choose here which editor should be used to display the result. If edi- 
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tor=none is set, the full result will not be displayed in an editor window at 
the end of the calculation. 

grouping If the togetherflag is set to 0, it could still be useful to first sum a 
few functions before integration. The number of bytes you set with group- 
ing= kbytes decides how many functions f*.f or f*.cc are first summed and 
only then integrated with the numerical integrator. If you set grouping=0 
all functions f*.f resp. f*.cc are integrated separately. In practice, a group- 
ing=0 has proven to lead to faster convergence and more accurate results. 
However, if you consider integrals which show large cancellations within the 
different functions f*.cc, it might be useful to use a groupingT^ 0. The log 
files *results*.log in the results directory contain the results from the in- 
dividual integration, where the user can see if there are large cancellations 
between the individual functions. 

language For one-scale diagrams or diagrams with purely Euclidean kine- 
matics language=fortran or language=Cpp (default) can be chosen, where 
the Cpp stands for C++. 

In all other cases, especially when using contourdef=True, language=Cpp is 
used, as the deformation of the contour which is needed for these problems 
is only implemented in C++. 

rescale If all invariants are very small or very large it is useful to rescale 
them to reach faster convergence during numerical integration. The rescaling 
(scaling out the largest invariant in the numerical integration part) can be 
switched on with rescale=l and switched off when set to 0. If switched on, 
it is not possible to set explicit values of any non-zero invariants in the 
Mathematica input file template* .m. 

contourdef For multi-scale problems resp. diagrams with non-Euclidean kine- 
matics, set contourdef=True (default is False). In this deformation 
of the integration contour in the form of Eq. ([6]) is done. In addition to 
the functions f*.cc to be integrated, auxiliary files (g*.cc) are written which 
serve to optimize the deformation for each integrand function. 

lambda Here, you can set the initial lambda A for the deformation of Eq. 
Without any knowledge about the characteristics of the integrand, lamhda=1.0 
should be a good choice. If the diagram contains mostly massless prop- 
agators and light-like legs, it can be useful to choose the initial A larger 
{e.g.lambda=5.0), in order to compensate for cases where the remainders 
of the IR subtraction lead to large cancellations for Xj — > 0. For diagrams 
with mostly massive propagators the initial lambda can be chosen smaller 
{e.g.la'mbda=0.1). 

smalldefs If the integrand is expected to be oscillatory and hence sensitive 
to small changes in the deformation parameter A, smalldefs should be set 
to 1 (default is 0). If switched on, the argument of each subsector function 
J-" is minimized. 

largedefs If the integrand is expected to have (integrable) endpoint singular- 
ities at Xj = or 1, the deformation should be large in order to move the 
contour away from the problematic region. If largedefs=l, the program tries 
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to enlarge the deformation at the endpoints. Further, the program performs 
a remapping to separate endpoint problems for Xj — >■ 1 from the ones for 
Xj — >■ if this flag is set to one. This increases the number of functions, but 
it can be very useful to improve convergence. The default is largedefs=0. 
In this case, the deformation parameter A is scaled down for each Feynman 
parameter. 

A. 2 Input for the definition of the integrand 

This Mathematica input file should be called *.m. The following parameters 
can be specified 

momlist If cutconstruct=0 is set in the input file, specify the names of the 
loop momenta here. 

proplist Specify the diagram topology here (mandatory). The syntax for cut- 
construct=l is described in Section 13. 4[ If cutconstruct=0 has been chosen, 
the propagators have to be given explicitly. An example propagator list 
could be proplist={k" 2-ms[l],(k+pl)" 2-ms[l]} with the loop momentum 
A;, the propagator mass m\ and external momentum pi. 

numerator If present, specify the numerator of the integrand here. If not 
given, a numerator={l} is assumed. Please note that the option cutcon- 
struct=l is not available in combination with numerator functions. 

powerlist As an option, the propagator powers (e.g. if different from one) 
can be set here. 

onshell Specify invariant replacements here. The kinematic invariants can be 
assigned values (e.g. ssp[l]— j-O) or relations between the invariants can be 
set (e.g. ssp[l]— >sp[l,3]). This option can not be used in combination with 
rescale=l. 

Dim Set the space-time dimension. The default is Dim=4-2 eps and the sym- 
bol for the regulator e must remain the same. 
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